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I. Introduction 



The presence of neutrino oscillations [1-2] has renewed interest in the question of 
oscillations within matter [3-5]. The early work ("MSW") of Refs. 3 and 4 solved the 
problem of propagation within a medium of constant density, and it is possible to treat 
nonconstant density by numerical means. It is nevertheless always useful to think about 
an analytic approach [6] in order to develop insight and understanding. Since the case we 
study here is that of a two-channel problem, some aspects of the methods we describe are 
also applicable to spins in varying magnetic fields and other two-channel order- 
dependent problems. 

For convenience we recall here the MSW results [3, 4], constant density We 
assume a two-channel approximation to neutrino mixing and give the amplitude T(t) for a 
neutrino beam of energy E passing through a medium of constant electron density N e 
given some initial neutrino flavor mixture y/(0), namely 

yKt) = T(t)ifKO). ' (1.1) 
(The time t can be interchangably viewed as the thickness x of the medium.) Here the 

\ V e 1 

two-vector weak state (i.e., flavor basis) is y/ (t ) = e . The transition amplitude T is 

T = cos^ / +/sin0 / cos(20')cr, -i sin 0' sin (26> ')(7 t (1.2) 



where the primed variables contain the effect of the matter: 
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A = 42G F N e . 

The mass parameter dm 2 = m 2 2 - mi 2 , where mf is the z' th mass eigenvalue, is positive. We 
recover the vacuum result, Cabibbo angle 9, for A = 0. 

This introduction sets some notation and recalls well-known results. Below we 
shall be concerned with variable density. 



II. Passage through a medium of variable density 



The Schrodinger equation for a two-family weak state y/(t) propagating through a 
medium of electron density N e is 

i*?- = (H F +W)yr, 
dt 



(2.1) 
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where, using the approximation m, « E and after pieces in the hamiltonian proportional 
to the unit matrix are removed — they lead only to overall phases, as we describe further 
below — we have 



Hi. = 



dm 2 



4E 



-cos 26 sin 26 
sin 26 cos 26 
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and 
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Here H F is related to the mass eigenstates by the Cabibbo matrix V = 



cos6> sin6> 
-sin6> cos6> 



H F = V ■ 



2E 



m 1 
ml 



■y 1 = vh d v\ 



(2.4) 



A solution to Eq. (2.1) will give us the amplitude T(t), through Eq. (1.1). 

It is convenient to approach the dynamical equation (2.1) by starting in the mass 
basis, using the transformed states (j>= Vty . We do so by multiplying Eq. (2.1) by V 1 , 
giving 

.dip 



i-L = (H D +V*WV)(l). 
Next we remove the factor H D by defining a new function % by 

Starting from Eq. (2.5) it is easy to see that this function obeys the equation 

i™L = e iH °'V*WVe- iH '>'£. 
dt 

In this expression, only the mass difference enters into the exponentials. To see this, 
define 



(2.5) 



. , in: +m\ 
E = E-\ — and A = 



AE 



dm 2 
~4E~ 
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In terms of these quantities, H D reads 



H D =E + (7 3 A 



The portion of Hu proportional to the unit matrix commutes through the quantities in Eq. 
(2.6) and cancels, leaving 



i d§_ = e iA< h tyt WVe -iA< h ,e 
dt 



(2.8) 



We next write out in 2 x 2 form the quantity multiplying ^ on the right side of Eq 
(2.8), using in particular the identity 





= cos a + 1<7. sin a 



e 



(2.9) 



We find 
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e^'V^WVe' 



' 1 + cos 20 e 2iA 'sin2e' 
e~ 2,A 'sin20 1- cos 20 
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1 ^ cos 26 e~ 

j- 2iA, sin2e -cos 29 
In the last step we have taken out the piece proportional to the unit matrix; again, it 
contributes only an overall phase to £ (This can be seen in a variety of ways. For 



- 2iA, sm2e 



(2.10) 



example, one can define a new function 77 = exp 



+ 



% and show that the 



equation for 77 is identical to that for % but with the last form in Eq. (2.10) on the right 
multiplying J].) We can finally transform away the term in Eq. (2.10) proportional to 03. 
We define the new function £by 



( i \ 
£=exp +—o 3 I (t)cos26 £ = 



exp 



-7(f)cos20 











exp 



7(?)cos20 



$,(2.11) 



where 



Ht) = \A{t')dt'. 



(2.12) 



The function £ obeys the simpler equation 
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A i(2Ar+/cos2e) • n, n 

Ae ' sin 20 



1 



(2.13) 



-Ae"' (2A ' +/cos2e) sin20 
.2 

The last line defines the matrix Pit), in terms of which we can write our formal solution 
to this equation. We remark in particular, for later use, that the 21 element P21 equals the 
complex conjugate of the 12 element Pn. We will below express the solution to Eq. 
(2.13) in terms of P n andPn*. 

The form that P(t) takes is easy to understand. The central matrix in the first line 

±-(T 3 /cos2e 

of Eq. (2.13) is a linear combination of 0\ and 02, while the external factors e 2 
take the form of a rotation about the 3-axis. Their effect is then to rotate the combination 
of Oi and 02 to give a different combination that lies at a different angle^ — that takes the 
form of the original combination but with the phase shifted. This is indeed what happens, 
as the explicit calculations that give us P(t) show. 
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Solution of the equation for £ It is clear from the matrix structure of Eq. (2.13) that the 
commutator [P(t), P(t')] ± 0, so that the solution to Eq. (2.13) involves an ordering. 
Formally, the solution is 

ft \ 



exp 



-ijp(s)ds 



C(0) S M(f)C(0), 



(2.14a) 



where the time ordering operator T ensures that the matrices P in the expression are 
ordered so that P with a later argument stands to the left of P with an earlier argument. 
We may write 

M{t) = f}^\ds]ds 2 --]ds n T[p{s l )P{s 2 )---P{s n )\, 
and it is a standard exercise to show that this is equivalent to 

M(t) = £(-i) n \p(s 1 )d Sl ]p(s 2 )ds 2 ---] 1 p(s n )d Sn (2.14b) 

1=0 o o 

Finally we can separate out the explicit matrix elements of M by using the matrix 
structure of P, namely P(t) = P x (t)Oi + P y (t)<72. This leads immediately to 



M n =\-\p n {A\p; 2 {f)dt'' 



dt'+. 



M n =-i\P n (t')dt' + i\P n (t') \Pn{t") \P n {f)dt" 



dt" 



dt'-. 



M 22 =M* U 



(2.15a) 



(2.15b) 



(2.15c) 



(2.15d) 

Alternatively, direct differentiation of the group of Eqs. (2.15) shows that it satisfies the 



necessary condition 



dM 
dt 



: -iPM 



The solution given by Eqs. (2.15) is essentially a power series in A. The order 
dependence of the result is contained in the fact that the integrals in M are nested. We can 
write the transition amplitude T in terms of M by undoing our series of transformations, 
leaving 



T(t) = V 



Af+-/(f)cos20 











;[ A/+-/(r)cos20 



M (t)V 1 



(2.16) 



Generally speaking the nested integrals in M are complicated, even for the case of 
constant density. (Of course numerical integration is always possible.) 

III. Second order equation 



Equation (2.13) represents two coupled first order differential equations. These equations 
can be rewritten as a single uncoupled second order equation, and the second order 
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equation lends itself to a variety of treatments beyond the formal solution we have 
already expressed. In terms of the explicit components of £, Eq. (2.13) reads 







"0 a 


"cr 




~aC 2 ~ 


_C 2 _ 




J8 0_ 


.c 2 _ 







where 



a = -iP\ 2 and jS = -iPi\ = -iP 



12 ' 



(3.1) 
(3.2) 



vc • 

One more derivative of, say, the upper component gives £j = a£ 2 + a£ 2 = — £\ + a/3£, , or 

a 



a 

Using the explicit expressions for a and fi, we find that 

- = — + i(2A+Acos20) and a/3 = 
cc -A 



(3.3) 



A sin 26 



\ 2 



Thus we have finally the second order equation 

— Ci+/(2A + Acos20)£ + 

i4 



1 > * 
— A sin 20 
2 



(3.4) 



It is worthwhile noting that the quantity /, Eq. (2.12), does not appear in this equation. 
We may take the required two boundary conditions to be £i(0) and £ 2 (0) = C, y (Qi)la(Qi) . 

1 dC,At) 

Once we find the solution for Ci(t), we have C 2 (t) = — — . We also remark here 

a(t) dt 

that we have verified that our formal solution to the equations for £ Eq. (2.15), satisfies 
Eq.(3.4). 



Recovery of constant density case. For the MSW case (constant A = Ao), reviewed in 
Section I, Eq. (3.4) takes the form 

Co d =0 

where 



b a = 2A + Ao cos 26 and c = ^ A,, sin 26 . 



If in addition we define 



then the solution to this equation is 



(0 = ^bl + 4c 2 , 



+ C, exp 



C 2 (0 = 5 2 exp |[«-^,] 



+ C 2 exp 



[-0) + b Q ] 
[-a -b ] 



where 



(3.5) 
(3.6) 

(3.7) 

(3.8a) 
(3.8b) 
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5, 



1 



ft) 



v 



ft) 



|c 1 (o)[«+M+c„c 2 (o) 



v 
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C 2 (o)[ft.+^ ]-c c 1 (o) 



v 

1 f 



ft) 



^(O^-^ + c^O) 
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This solution allows us to identify M(t) through Eq. (2.14) and hence the 
transition amplitude T through Eq. (2.16) adapted to constant A, namely 



T(t) = V 



-(2A+/l cos2e) 







|(2A+/l cos2e) 









MV f =V 


e 2 










e z 



MV 1 



(3.9) 



T = — V 
2(0 



■imt/2 



e 

The result of the exercise is 

~(ft) - b ) e im " 2 + (ft) + b ) e- ia " 2 -2c (e ia " 2 - e^ 11 ) 

-2c (e imn - e- im ' 2 ) (ft) + b ) e to,/2 + (ft) - & ) e 

Then specific calculation shows, for example, 

T n (= T eM ) = -± A sin 26 sin (ft)? / 2) . 

This can be put into the canonical form of Ref. [4], -z'(sin<p,„)(sin26> m ), where the subscript 

dm 2 ff 

m indicates the propagation is in material, and where q) m = — t , if we identify 

4/: 



V 1 (3.10) 



(3.11) 



sin 20, = ■ 



2|A|sin20 



sin 29 



sin 29 



ft) 



I 4EA ^2£A^ 
1--— r cos 20 + 



dm 2 



y (5m 2 y 



and 




(3.12a) 



8m 2 ff = 2E(0 = dm 2 Jl — 



4EA 
(5m 2 



cos 20 + 



2£A 
(5m 2 



(5m 2 



cos 20 - 



2£A ^ 
dm 



+ sin 2 20 .(3.12b) 



Indeed, in terms of these new variables the full transition amplitude is 

T = cos <p m + (i sin <p m cos 20 m ) o z - (i sin (p m sin 20 m ) o x . (3.13) 

This is the full canonical form described in Section I for propagation in a medium of 
constant density. 

Adiabatic Expansion. If the factor in Eq. (3.4) that contains the derivative of A is small 
compared to the other factors, one can make a systematic adiabatic expansion [6zz] in 
terms of it about the th order (MSW) answer. To do so, it is useful to recast the solution 
technique somewhat. We shall first take the starting point of the neutrino beam, at t = 0, 
to specify the constant background level of the material density factor, i.e., A(0) = A . We 
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leave the boundary conditions £i(0) and ^(0) unspecified for the moment but remark that 
using Eq. (2.13) the boundary condition for £ 2 (0) gives us alternatively a condition for the 

derivative of £i at t = 0:^(0) = ~ iP n (0)£ 2 (°) = ~ A> ( sin 20 )C 2 (0) • 

We see from our earlier solution of the constant density case (Eq. 3.8) that the t- 
dependence is contained in a pair of phases. An alternative way to derive these phases in 
the constant density case is through a solution ansatz of the schematic form 

Ci = *oexp(iSo(0), (3.14) 
where R is constant and where So(t = 0) = 0. The real and imaginary parts of Eq. (3.5) 
lead to the following equations for So(t): 



dt 



A L _ 2 

at 



d 2 S, 



(3.15) 



^ = 0. 



dt 2 

Together with the vanishing of So at t = 0, these imply that 

S (t) = Xt (3.16) 

where 

X 2 -b^-cl=Q (3.17) 
The solution of the quadratic equation gives 

X ± =^(b Q ±co), where co = ^jb 2 + 4c 2 . (3.18) 

Thus, as indeed Eq. (3.8) shows, a better ansatz for the solution is 

C 1 {t) = R + e a+t + R_e aj (3.19) 
The boundary conditions for give us immediately R+ + R-= £i(0) and 

R + X + +R_X_ =--^j4 £ 2 (0) sin 20 ; in turn these last two relations determine 

r _ C 1 (0) _ ^C 1 (0) + 4,C 2 (0)sin2g 

± 2 + 2co 
Having reviewed the th order (constant density) problem, we go on to include 
time (distance) dependence in the material density. We accordingly write the input 
density as 

A(0=A (l+/i(0), (3.21) 
where f\{t) « 1 for all t in the problem and/i(? < 0) = 0. We extend our ansatz for the 
solution to the form 

C, (0 = R + (0 f3 + e iS + (t) +R_ (t) f3_e iS ( '\ (3.22a) 
where R ± (t) = R ± (l + pf (tj) (3.22b) 

S ± (t) = X ± t + af(t) (3.22c) 

The quantities with subscript "1" are all small; moreover, pi -(0) = of (0) = 0. We also set 

b(t) = b + bi(t), where bi(t) = A fi(t) cos20 (3.23a) 
c(t) = co + ci(t), where a(t) = (A fi(t) sin20)/2. (3.23b) 
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We now insert our ansatz into Eq. (3.5). The spirit of the adiabatic expansion is to keep 
only first order terms in quantities with the subscript "1." In addition, we insist that the 
coefficients of exp(zS+(f)) vanish separately. 

The real and imaginary parts of the coefficients of exp(z'5±) give respectively 

d 2 pf d<77 „ , _ . . 

— ^- + co—±- +K + f 1 = (3.24a) 
at at 

^±0)^-^ = 0, (3.24b) 
dt dt dt 

where 

K± = A o (A±cos20 + c o sin20) = 2A±(-2A ± (0)12. (3.25) 

Equations (3.24) contain only derivatives of the functions we seek, so they are in 
fact two coupled first order equations for the functions 

± do* ± dpf 

v l = — - and u l =— —. (3.26) 
dt dt 

To the equations for vr and u\~ we add boundary conditions that follow from dC,\ldt = 0, 
namely vr(0) = = wr(0). As we shall see, these boundary conditions guarantee that 0\ 
and pi remain small (i.e., 0(f\)). 

We decouple the two equations for vi* and wi* by taking one more derivative of, 
say, Eq. (3.24a), giving a single second order equation for u*: 

^t + G) 2 uf =2AK~, (3.27) 
dt dt 

where for the coefficient of df\ldt we have used ±0)X± - K± = 2A. This equation is solved 

by a standard Green's function G(t - t') that satisfies 

d 2 G(t-t') , . ,i - / /x 

K - — J - + (d 2 G{t-t) = 8(t-t) 

dt 

with, as causality suggests, G(x) = for x < 0. The Green's function required is 

G(x) = — 0(t)smeot. (3.28) 

ft) 

In terms of this function we have 



uf (t) = a ± s sin at + a* cos cot + 2AA ± I dt'G (t - 1') — ^ - 



= a: sin ft)? + a* cos a)? + 2AA + — f cfr' — - - sin (ft) (t - 1')) 

~ col dt 

(The lower limit reflects the fact that/i is identically zero for negative argument.) Once 
we have u-^ we can get vi* from Eq. (3.24a): 



(3.29) 



CO 



= + 



(3.30) 

^ + K ^ 

cos ft)? - a~ sin ft)? H — - /, (?) 

ft) 
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At this point we can get a s ± and from the boundary conditions for vi^O) and wi^O), 
which determine a s ~ = a c ~ = 0. (It is worth noting that since these quantities are not 
proportional to/i, the only way for them to be small is to be zero.) In turn, we have 
finally 

^ li'sinM*' ' 
CO 



p (t) = 2AA j jsin ( ( o{i'-i"))'' 







df 



(3.31) 



(7 1 ± =±-^ + f/ 1 (0* / 
CO J n 



(3.32) 



The phase function S±(f) takes on a suggetive form if we use the identity K± = ±CoA+ - 
A(bo ± CO), in which case we can write 



s ± (0 = a, J (i + /x (OK t — J /i (0 ■ 

W 



(3.33) 



The first term integrates the material density. 



An example: We take a linear variation, f\ = qt, together with the condition that 
the beam is pure at t = (which translates into £i(0) = -sin6> and £ 2 (0) = cos6>). Then 
we have immediately 



S ± (t) = X ± 



t + -qt 
2 



+ - 



2AA ± 1 
co 2 



qt 



and 



R ± {t) = R ± 



l + q 



2AA ± 



co 



f 1 . 
? sin CO? 

CO 



(3.34a) 



(3.34b) 



In Figs. 1 and 2 we plot some probabilities associated with the resulting amplitude for 
some representative values of the parameters. 
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Figure Captions 

Figure 1 . Probability, as calculated in the adiabatic approximation described in the text, 
of V e as a function of time from production as a pure Vu at t = 0, in a medium with density 
factor Ao(l + qt), where Ao = 6 x 10 9 cm -1 = 10~ 13 eV corresponds to Earth-like density. 
We assume the primary mixing angle is = 0.7 and that the difference of the square of 
the neutrino masses is 5 x 10~ 6 eV 2 . The factor A = -7.9 x 10~ 14 eV, a value for which the 
neutrino energy lies around the MSW resonance value, an energy of roughly 20 MeV. 
The horizontal exis is in units of 10 14 eV -1 ; note that for q = -2 x 10~ 15 eV, qt = -0.2 at t 
= 10 14 eV" 1 . 

Figure 2. Same as Fig. 1, but with the factor A = -7.5 x 10~ 15 eV, a value for which the 
neutrino energy lies roughly ten times higher than that corresponding to the MSW 
resonance value. 
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